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ABSTRACT 

The theory and procedures for determining the characteristics of pressure 
oscillations in rocket engines with prescribed burning rate oscillations are 
presented. Pressure and velocity oscillations calculated using this procedure 
are presented for the SSME to show the influence of baffles and absorbers on 
the burning rate oscillations required to achieve neutral stability. Results 
of calculations to determine local combustion responses using detailed 
physical models for injection, atomization, and vaporization with gas phase 
oscillations in baffled and unbaffled SSME combustors are presented. The 
contributions of the various physical phenomena occurring in a combustor to 
oscillations in combustion response have been determined. 

INTRODUCTION 

Any review of liquid rocket combustion instability model development in the 
United States prior to 1972, inevitably leads to NASA SP-194 This 
publication is a compendium of the combustion instability modeling efforts by 
industrial, academic, and government teams from across the nation. Their 
efforts were aided by the considerable financial support these groups received 
during the 1960 's due to the perceived threat of combustion instabilities to 
achieving national space goals. The work of these groups also benefitted from 
the considerable experimental studies of combustion instability phenomenon 
including liquid rocket combustion tests. 

After 1972, financial support for combustion instability research was 
drastically reduced and very little effort was devoted to the subject. In the 
interim, major advances were made in the computational capability of computers 
and in numerical solution techniques for non-linear partial differential 
equations. This increase in capability has permitted investigators to make 
calculations that were only dreams in the early 1970's. Increased 
computational capabilities also allowed the development of new quantitative 
spray measurement techniques . 

Atomization was shown to be an important mechanism in combustion instability 
in numerous experimental studies. Experimental studies 2-4 clearly illustrate 
the drastic spatial and temporal changes that occur in the atomization process 
during an instability. To this day, the quantification and simulation of the 
atomization process remains a considerable challenge to the understanding and 
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control of combustion instability in liquid rocket engines. 

STATUS AS OF 1972 

The majority of the modeling efforts prior to 1972 used linear or nonlinear 
acoustics but very cursory combustion models. Modeling efforts representative 
of this era were made by Crocco and his colleagues at Princeton 5-8 , Culick 9 , 
and Priem 10 . These models were used to gain an understanding of instabilities 
and guide engine design changes intended to achieve stability. In general, the 
differences between the chamber response obtained with linear and nonlinear 
acoustic models were small 6,7 . Some models for chamber acoustics included the 
separate effects of baffles (radial baffles or hub baffles but not 
combinations of hub and radial baffles) and acoustic absorbers but not the 
combined effects of baffles and absorbers. With these models of chamber 
acoustics, designers had a greater ability to predict the need for and proper 
design of acoustic absorbers. 

The combustion models used in the frequency domain codes related the local 
instantaneous combustion rate to the local pressure using two constants (i.e. 
n and tau or a complex combustion response) . Reardon 11 did correlate time lag 
parameters with engine design and operating parameters. These correlations 
were also used as a guide to obtaining stable engine configurations. 

Individual combustion processes such as injection, atomization, and 
vaporization were too computationally intensive to model and in the case of 
atomization detailed models were not available. Movies of the combustion zone 
in experimental combustors 2 have shown the highly nonlinear nature of the 
combustion process during an instability. During a full-blown instability, it 
almost appeared as if a detonation wave was traversing the combustion zone. It 
is not suprising that the simple two constant combustion models could not 
provide a detailed understanding of instability or a quantitative evaluation 
of the influence of design and operating changes on stability. 

Numerical studies 12,13 were conducted to examine the response of individual 
combustion processes (injection, atomization, vaporization, and chemical 
reaction) to idealized chamber acoustic waves. These acoustic waves 
corresponded to the resonant frequency of cylindrical chambers with no 
stability aids, axial velocity profiles, and ideal nozzles. The studies 
demonstrated that all of the individual processes were capable of driving an 
instability. The result of these studies were correlated to provide n-tau 
combustion models. None of the frequency domain models of this era attempted a 
comprehensive model that tracked the propellant from the injector until it was 
completely burned. 

Simple one and two dimensional time domain combustion instability models 14 ” 16 
were developed prior to 1972. Although these codes typically modeled 
injection, a mechanistic treatment of atomization was omitted. The combustion 
rate was usually modelled as an empirical function of pressure. None of these 
early time domain codes had a "complete” combustion model. The codes were able 
to demonstrate the existence of a minimum disturbance amplitude necessary to 
initiate an instability. The wave characteristics calculated with these codes 
were similar to those observed from dynamic stability tests. These codes had 
considerable problems with numerical stability, particularly the two 
dimensional codes . 
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DESIGN TOOLS FOR THE 1990 's 


The goal of current stability code development efforts is to produce a 
collection of analysis tools to allow a designer to assess the stability 
margin of candidate combustor designs. This collection of codes may include 
fast, acoustic n-tau codes for preliminary design screening. Fourier time 
expansion and simple time domain codes could be utilized for stability margin 
evaluation and stability aid design for candidate combustor designs. Rigorous 
time domain codes could be used for final design evaluation. The features of 
these codes are summarized in Table I. and are discussed below. 

Any analysis of combustion instability in rocket engines is only as accurate 
as the specification of the boundary conditions for the problem. 

Unfortunately, the boundary conditions for stability analyses in liquid rocket 
engines are highly nonlinear and not very well understood. The rate and 
location at which injected propellant is atomized and the resulting size and 
velocity of the droplets is crucial input to any instability analysis. The 
subsequent behavior of the propellant as it burns or undergoes further 
breakups is critically linked to the atomization boundary condition. This 
boundary condition is also part of a highly nonlinear feedback loop in which 
the combustor gas field influences atomization and atomization influences 
burning and the combustor velocity profile. Accurately modeling this phenomena 
is the paramount challenge of rocket combustor modelling. 

Typically, a shear coaxial element is modeled using correlations developed 
from cold flow tests. These correlations provide a relationship between global 
parameters (gas velocity, liquid flowrate, jet diameter, etc.) and drop size 
or atomization rate. Unfortunately in many models, these global correlations 
have been applied in a subgrid fashion. The overall significance of droplet 
coalescence- and secondary breakup phenomena has not been demonstrated. 
Secondary breakup is typically modelled using a spring mass oscillator analogy 
or linear surface wave growth rates even though the majority of experimental 
evidence suggests a violent stripping phenomena. One can imagine the 
difficulty in trying to model the drop size produced by an impinging element 
subject to a high amplitude transverse wave. 

The biggest advantage of time domain codes is that they can provide a complete 
time history of any parameter of interest. These time histories allow for a 
direct comparison with engine test results. Physical submodel development and 
code validation is also facilitated by use of these time histories. 

Interaction between different modes of oscillations is handled automatically, 
within the filtering limits of the finite difference method. 

Updating of one and two dimensional time domain codes to take advantage of the 
increase in knowledge of the stability characteristics of various finite 
differencing techniques is underway. However, most instabilities due to 
injector design or the wave dynamics are inherently three dimensional. A fully 
three dimensional simulation of the SSME (manifolds, baffles, absorbers, and 
600 injection elements) is beyond current computational capabilities, even if 
atomization could be modeled adequately. Computational time remains a major 
drawback of the time domain approach. 
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TABLE I 


THREE APPROACHES TO INSTABILITITY MODELING 


Acoustics - 

Simple Combustion 

FEATURES : 

-linear amplitude or 

mode series expansion 

-continuity and momentum 
conserved 

-simple combustion 
two parameters 

-stability aids modeled 

ADVANTAGES : 

- identifies problem 
modes 

-results obtained 
quickly 


Time Domain 
CFD 

-time history all 
parameters . 

- continuity , momentum , 
and energy conserved 

-mechanistic combustion 
model 

-stability aids modeled 


-includes all nonlinear 
effects 

-direct comparison to 
stability tests 


Fourier Time 
Expansion 

-Fourier expansion 

in time 

-continuity and 
momentum conserved 

-mechanistic 
combustion model 

-stability aids 
modeled 


-determines 

margin 

-simulates resulting 
oscillation 


DISADVANTAGES : 
-poor modeling of 
combustion 


-computationally intense -no real time 

history generated 


-difficult to determine 
instability source 


-different runs for each -no interaction 
disturbance level, type between modes 


The Fourier time expansion (FTE) approach is something of a hybrid between the 
time domain, CFD approach and the analytical, acoustic modeling approach. In 
the FTE approach, the thermodynamic variables and velocity components are 
expanded in a Fourier time series and substituted in the appropriate 
conservation equations. The resulting equations are differenced and solved 
numerically to determine the gas phase oscillations. All the basic combustion 
processes are modeled as they would be in a time domain approach and then 
converted into complex responses for use with the FTE approach. Solutions are 
obtained at different oscillation amplitude levels for the complex frequency 
(observed frequency and decay rate) . The oscillation amplitude at which the 
decay rate changes from negative to positive is the minimum disturbance 
amplitude required to excite that mode of instability. The oscillation 
amplitude at which the decay rate changes back from positive to negative is 
the equilibrium or "limit cycle" amplitude of the instability. A disadvantage 
of the FTE approach is that it does not permit having oscillations present 
with different basic frequencies, i.e. chug and a high frequency oscillation 
or two transverse modes with different frequencies that beat. Another 
disadvantage is that no time history of flow parameters is calculated directly 
for comparison with experiments. 
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THEORY FOR THE FOURIER TIME EXPANSION 


In the remaining portion of the paper, the theory for the Fourier time 
expansion (FTE) model will be discussed and calculations made using the method 
will be presented for an SSME combustor. The results of the FTE approach was 
developed for incorporation into HICCIP (High Frequency Injection Coupled 
Combustion Instability Program) 24-26 . The code described in ref. 24 was 
restricted to using wave characteristics derived from a pseudo 3D theory. The 
code assumed uniform combustion response at all radial and tangential 
positions for a given axial position. Also, the wave model used could not 
account for the presence of baffles and discrete absorbers. Therefore, a 3D 
model for the chamber gas oscillations was required for HICCIP to predict the 
characteristics in a typical rocket engine. 25 ' 26 

The gas phase oscillations calculated with the 3D FTE approach are used to 
determine burning rate oscillations throughout the combustor. This is an "open 
loop 11 calculation where burning rate oscillations are assumed and used to 
calculate gas phase oscillations. The gas phase oscillations are then used to 
calculate burning rate oscillations . If the final calculated burning rate 
oscillations agreed with the assumed burning rate, a consistent, complete 
solution would have been obtained. To close the loop, and obtain the complete 
solution would require iterating the procedure described above. Work is 
continuing toward completing the "closed loop" model. 

The basic conservation equations used to define flow in a cylindrical chamber 
are : 

CONSERVATION OF MASS: 




CONSERVATION OF MOMENTUM: 


3( P V) +V- P VV +S rVP= Wi lq 


ADIABATIC IDEAL GAS: 



Where : 


P 

V 

P 

t 

W 


g 

V 


gas density, lb/in 3 
gas velocity vector, in/s 
pressure, lb/in 2 
time, s 

local burning rate , lb/in 3 s 
velocity of the injected liquid, in/s 
conversion constant, lbm in/lbf s 2 (386.4) 
space derivative, in -1 
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It is assumed that viscous and gravitational forces are negligible. To 
simplify the analysis, it is further assumed that the radial and azimuthal 
steady or average flows are minimal and are therefore zero. This avoids 
requiring a 3D steady state or time average solution to determine radial and 
azimuthal flow and pressure profiles. 

As in HICCIP , all the variables (P, p, V z , V r , V a , W) are expanded 
in a Fourier time series as given by: 


F=F 0 +F 1 e iwt +F 2 e 2iwt + +F D e niwt +. . . 


Where: F 0 is the time averaged variable 

F n is a variable in all three dimensions 

Consistent with the above assumptions: 

p 0 & P Q vary only in the axial direction 

V ro 6c V ao are zero (no radial or tangential average flow) 

V zo & W Q are variable in all three space dimensions. 

All of the equations are nonlinear. In the remaining portion of this paper, 
only the terms through F x are used as these terms are not influenced by the 
F 2 , and higher n terms (since low frequency cross products produce higher 
harmonic frequencies) . After determining the F 1 values the equations could 
again be solved using the known F 0 and F x values to find the F 2 values. This 
process could be continued to determine as many F n values as required to 
"completely and rigorously" define the oscillations. Solutions for the F 1 
values will define the complex frequency of the oscillations (both the 
observed frequency and the damping or growth rate) . 


Substituting the Fourier time series into the conservation equations and 
retaining only the terms containing "F 0 and F x " yields the following set of 
equations : 

CONTINUITY: 


1 d(r 9o V Zi ) x d(p Q V ) d( Po V ) d( 9l V ) 

10 )p i + T ai — '~ ae * d z * Si " 




AXIAL MOMENTUM: 


Eq.l 


w . .2 atr -W ♦ airp f ■ v - ) ^ 

ko zi n z Q dz r dr r dd 


+ Y- 


3pi 


Eq. 2 


RADIAL MOMENTUM: 


d(9 0 V z V Zi ) Po 


ioPoVrS dF~^~ +y 


dp i 
Po dr 


=0 


Eq. 3 


6 


AZIMUTHAL MOMENTUM: 



Eq.A 


ADIABATIC IDEAL GAS: 



Eq. 5 


The boundary conditions for the combustor are: 


Chamber Wall 


V Xi = y RspWall P x 


where "RspWall" is zero for hard walls and is an input value for wall absorbers 
(RspWall is ratio of the radial velocity to local pressure) . 

Center of Chamber ( r = 0 ) 

All radial fluxes are zero as the area is zero. 

Injector Face 


where "Rsplnj " is zero for a solid injector face and is input for absorbers on 
the injector face. 

End of Combustion Chamber 


where "RspNoz" is the nozzle response obtained from another program, theory etc. 
as desired. For this paper it is assumed to be zero (short distributed ideal 
nozzles) . 

Baffle surfaces 

Radial Baffle V al = 0 Baffle is at an azimuthal momentum cell face. 

Hub Baffle V rl = 0 Baffle is at a radial momentum cell face. 

Equations 1-4 result in a very sparse matrix when using finite difference 
equations as illustrated above. Using upwind differencing for fluxes due to 
steady flow the number of terms in each equation are: 


V Zi = y Rsplnj p 1 


V Xx = Y RspNoz f> 1 
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continuity 
axial momentum 
radial momentum 
azimuthal momentum - 


8 terms 
13 terms 
4 terms 
4 terms 


The total number of equations and variables to be solved is determined by the 
number of cells used in each of the three dimensions as given by: 

equations = variables — 4 * Ni * Nj * Nk 

Where: Ni number of radial cells 


As an example, for Ni=10, Nj=12, and Nk=20, which adequately defines an 
unbaffled chamber with the fundamental transverse mode there are 12,000 
equations with complex variables. Each equation has an average of 7 unknowns. 
The "W" source term is specified for each cell and is determined using burning 
models as in HICCIP. 

One dimensional or coarse two dimensional solutions can be obtained with a 
banded Gaussian elimination procedure. For problems with larger grids , the 
resulting sparse matrix equations are solved using a conjugate gradient 
algorithm modified for complex variables. Since the matrix is neither symmetric 
nor positive definite, a quadratic form of the error is minimized. This 
quadratic form effectively has the square of the condition number of the 
original matrix. As a result, it is necessary to use double precision on 
machines with 32 bit word sizes. The current algorithm does not incorporate 
preconditioning. Crude attempts to implement a "multigrid" solution procedure 
did not reduce computation time. 

ID and 2D MODEL 

It is also possible to model the oscillations as ID or 2D phenomena. This is 
accomplished by assuming that the oscillations in the radial and/or azimuthal 
directions can be specified by an analytic wave equation solution in the desired 
direction. With this assumption, the conservation equations in that direction 
are not necessarily satisfied. The wave equation is: 


Nj number of azimuthal cells 
Nk number of axial cells. 


(z) e iwC e i7J0 J r r (/nr) 




and the variables are determined by: 
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For one cell, in the radial direction the radial derivatives in equations 1 to 5 
are given as follows: 


l J_ 
r dr 


(r Po v ri ) lJL( z j a lmz))/-£ ) ]v ri 


T fz 


(rJ n (mr) /£j a (mr) lv Xi 


P 0 9p t _ YP 0 
Po ^ Po 


1-^: ( rJ D (mz ) ) /J n (mr) lp x 


where the Bessel derivative terms are determined from Bessel function tables (at 
a radius equal to half the chamber radius) for different Bessel Function indices 
M m and n" which describe the mode in the radial direction and the tangential 
mode (n) . The calculations then correspond to a very thin annular ring at half 
the distance to the chamber wall. 

Similarly for only one cell in the azimuthal direction the azimuthal derivative 
in equations 1 to 5 are given as follows: 


1 dp„V ai 

i de 



i a<PoVV 
r de 


Po*'*, 


inV a 


1 d(PoPi) 

i d6 



where n specifies the number of modes in the azimuthal direction. The 
calculations then correspond to a very thin "pie" section cut out of the 
chamber. 


EXERCISING THE MODEL ON THE SSME 


To demonstrate the capability of the Fourier time expansion (FTE) approach, 
calculations were made for the SSME at nominal operating conditions with and 
without baffles and absorbers for different oscillation modes. The configuration 
used to represent the SSME is presented below: 


SSME CONFIGURATION 


Chamber Diameter 

17.74 

in. 

Cylindrical Length 

15.079 

in. 

Hub Baffle Length 

2.0 

in. 

Hub Baffle Diameter 

8.87 

in. 

Radial Baffles 

5 

Blades 

Radial Baffle Length 

2.0 

in. 

Chamber Pressure 

2844 

psia 

Oxidizer Flow 

834.0 

lb s/s 

Fuel Flow 

225.3 

lbs/s 

Absorber Wall Gap 

0.36 

in. 

Absorber Tuned Freq 

3T 

Mode 

Speed of Sound 

5183 

ft/s 


SSME GAS PHASE OSCILLATIONS 

A time averaged combustion profile was used similar to that obtained with a 
vaporization limited model (50% combustion completed 2.3 inches from the 
injector, 90 % in 7.5 inches, and 99 % completion within the chamber) . The time 
averaged combustion profile remained constant. This is not a realistic 
assumption as the gas oscillations do influence the time averaged burning rate. 
The influence of these oscillations will be evaluated when the FTE approach is 
combined with the combustion models in HICCIP. 

The definition for combustion response is: 

Wi / W G 

CombRsp = 

Pl / Po 


where : 

W x is the perturbation in burning rate, lb/(s in 3 ) 

W Q is the local average burning rate, lb/(s in 3 ) 

P : is the perturbation in pressure, psi 

P G is the local average pressure, psi 

The combustion response is a vector in the complex plane. The real portion of 
the combustion response is the magnitude of the nondimens ional burning rate 
oscillation that is in phase with the pressure. The imaginary portion of the 
combustion response indicates the phasing of combustion relative to pressure. A 
positive imaginary response indicates that the combustion is leading the time 
dependent pressure, a negative response indicates combustion is lagging. 

Burning rate oscillations were obtained by using a combustion response that did 
not vary spatially. The simplif ications of the burning process, while not 
rigorous, allow a demonstration of the FTE approach and indicate how some of the 
combustor design variables influence stability in the SSME. 
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The combustion response needed to drive various waves in the SSME were 
calculated iteratively. For a given frequency, a combustion response was 
assumed and a burning rate oscillation was approximated using an analytical 
acoustic wave solution. This initial guess at the burning rate oscillation was 
used to determine oscillating pressure and velocity profiles in the combustor. 
With the calculated pressure profile and assumed burning rate oscillations, a 
mass averaged combustion response was calculated. The mass averaged combustion 
response was then used to calculate a new burning rate oscillation. This process 
was repeated until the difference in the burning rate was less than 0.003 
between successive calculations. 


With the assumptions and procedure outlined above, a combustion response for any 
frequency can be determined. The combustion response was calculated (using a 
grid size of 4 radial, 10 azimuthal, and 8 axial cells) for the IT mode and the 
results are shown in Fig 1. The complex combustion response is plotted as a 
function of frequency between 1860 and 2060 Hz. Three curves are shown in Fig. 1 
corresponding to damping rates of +628 s” 1 , 0 s" 1 , -628 s" 1 . As the frequency is 
increased, the combustion response changes from combustion lagging pressure to 
leading pressure and crosses the in phase (1.39 real and zero imaginary) plane 
at 1960 Hz with neutral damping. Normally, this is considered the "natural 
frequency" for the IT mode in this engine. However, depending on the combustion 
response, oscillations in the engine can occur at other frequencies as shown in 
Fig 1. Larger combustion responses make the engine unstable, lower values result 
in a stable engine. Negative imaginary response lowers the tuned frequency from 
the natural and positive imaginary responses increase the tuned frequency. 
Similar plots can be obtained for all the modes possible in the chamber. The 
mode type is determined by the initial assumption for the pressure profile or 
burning rate oscillations in the chamber. 

To illustrate how the FTE approach can be used to determine the influence of 
design parameters on stability and stability modes, calculations were made to 
determine the frequency and combustion response for "neutral stability at the 
natural frequency", i.e. where the combustion response has no imaginary value 
and the oscillations have zero growth or decay rate. The results of these 
calculations for different stability aid configurations and modes are presented 
in Table II. 
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TABLE II 


NEUTRAL STABILITY CONDITIONS FOR SSME 


Mode 


Frequency 

Response 


UNBAFFLED 

AND 

NO ABSORBER 

(Grid of 

8*1*16) 


1 L 


1966 


1.69 


1 T 


2057 


1.69 


2 T 


3407 


1.65 


3 T 


4686 


1.64 


UNBAFFLED 

AND 

ORIGINAL SSME 

ABSORBER 

(Grid of 

8*1*16) 

1 L 


1955 


1.73 


1 T 


2045 


1.63 


2 T 


3366 


1.50 


3 T 


4680 


2.63 


UNBAFFLED 

AND 

TUNED TO SPECIFIC FREQUENCY (Grid 

of 8*1*16) 

1 L 


1966 


1.94 


1 T 


2052 


1.90 


2 T 


3403 


2.51 


3 T 


4680 


2.63 


BAFFLED AND ORIGINAL SSME ABSORBER 

(Grid of 8*20*16) 

1 L 


1955 


1.73 


1 T 


1957 


1.35 


2 T 


3023 


1.08 


3 T 


4305 


2.26 



Table II illustrates how the acoustic absorber improves stability in the SSME. 
The original absorber design is tuned to the 3T mode and increases the 
combustion response of the 3T mode for neutral stability from 1.70 to 2.63 which 
is equivalent to going from neutral stability to a damping rate of -628 1/sec in 
Figure 1. This damping rate would decrease the amplitude of a 100% chamber 
pressure pk-pk wave to the allowable 10% amplitude in 3 . 6 milliseconds. The 
damping rate as required by ref. 28 is 18 milliseconds . Therefore, the SSME 
absorber design provides damping 5 times that required for the 3T mode. 


The original absorber design provides a very small change in frequency and 
response for the IT and 1L modes because the absorber is tuned to a much higher 
frequency. The 2T mode is actually made more unstable by the absorber because 
it reduces the frequency and decreases the chamber response. Re tuning the SSME 
absorber to the IT and 1L modes by decreasing the absorber opening in the 
chamber from 0.36 in. to 0.12 in. and increasing the aperture length from 0.5 to 
1.0 inch will produce a large increase in stability, more than enough to meet 
the damping requirements of ref. 28. The increase in stability at the IT and 1L 
frequencies is less in comparison to the stability of the original design at the 
3T because the absorber open area has been reduced to obtain tuning at these 
lower frequencies. 

Adding the baffle to the SSME changes the tuned frequencies and decreases the 
response required for neutral stability. The frequency change is 105 Hz for the 
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IT mode and 380 Hz for the 2T mode. The reduction in response is also greater 
for the 2T mode. 

These calculations demonstrate that oscillations can occur in baffled chambers. 
Oscillations with either standing or spinning wave characteristics are possible. 
With the baffle the wave is distorted within and near the baffles, but an 
acoustic like wave does exist. Downstream of the baffles the wave is similar to 
the acoustic wave for a pure cylinder. Near the baffles they become distorted 
and within the baffle compartment the oscillations are predominantly axial with 
small variations in the radial and azimuthal directions. The decrease in 
response occurs because of the wave distortion produced by the baffles coupled 
with the steady flow profile within the chamber. The baffles reduce the loss 
associated with the steady flow convecting spinning oscillations out of the 
chamber by producing axial oscillations . 

To test the capability of the FTE approach to predict observed acoustic 
phenomena, calculations were performed to compare with the results published by 
Wieber 29 and the Gems ip Program 30 . Good agreement with this data was obtained. 
Calculated results for the SSME were somewhat dependent on grid size, especially 
for the baffled chambers . 25,26 


SSME BURNING RATE OSCILLATIONS 

Local burning rate oscillations resulting from the gas phase oscillations 
described above were calculated using the HICCIP model described in ref 24. 
These calculations were performed at the resonant frequency (no imaginary 
combustion response and no damping rate) for the SSME chamber and for the SSME 
with stability aids removed. The HICCIP model includes: 

1. Fuel flow oscillations 

2. Oxidizer flow oscillations 

3. Atomization in the recess and chamber using Faeth (ref 17) 

4. Vaporization using an onion skin model (ref 20). 

With these models, the propellant is followed from the time it enters the fuel 
or oxidizer inlet dome, through the injection element, recess, and combustion 
chamber (as a function of time and space) . Mass is always conserved throughout 
the combustion process. The processes are dependent on the local conditions as 
required by the available equations. 

In order to make the burning rate oscillation calculations, a proper selection 
of the constants in the atomization models for the rate of atomization and drop 
size is required. The atomization rate constant was selected so that 50% of the 
oxidizer was atomized in the recess. This resulted in an atomization plane (the 
effective location where vaporization starts) at the injector face. The drop 
size constant was selected to produce a steady state gas velocity profile 
similar to that used in the gas phase oscillation calculations described above 
(99% vaporized in chamber). This resulted in a steady state drop size of 507 
microns . 

The oscillations in injector flow rates as given by a flow response for 
different rows in the injector was: 
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TABLE III 

INJECTOR FLOW RESPONSES 


Row Fuel Response Oxidizer Response 


1 

-5.9 

+ 1.2 * i 

-0.42 

- 

0.30 * 

i 

2 

II 

It 

-0.53 

- 

0.40 

ti 

3 

It 

ti 

-0.70 

- 

0.46 

it 

4 

11 

ti 

-0.93 

- 

0.44 

it 

5 

It 

if 

-1.16 

- 

0.25 

ti 

6 

n 

it 

-1.22 

+ 

0.08 

ti 

7 

it 

n 

-1.06 

+ 

0.34 

tt 

8 

it 

ti 

-0.83 

+ 

0.44 

it 


The large fuel response is due to the very low fuel injection pressure drop (305 
psi) . The variation in oxidizer response is due to the different lengths of the 
oxidizer tubes connecting the oxidizer manifold to the combustion chamber in the 
SSME injector. Oxidizer tube length varied from 5 in. near the center to 9 in. 
at the chamber wall. 

Local combustion responses are shown in Figure 2 as a function of axial position 
for five different radial and tangential positions. In order to compare the 
local combustion rate oscillation at different locations they are 
nondimensionalized as follows 

n- P »L y » P > 


where 

is the local oscillating burning rate 
W ol is the average burning rate in the first axial cell 

P Q is the average chamber pressure 

P n is the oscillating pressure in the first axial cell 

Thus cells with very low burning rates and/or pressure amplitude do not 
contribute to the oscillations and are reduced proportionally in combustion 
response. The responses are referred to as "normalized local combustion 
responses" in Figure 2. This weighting eliminates large values of the combustion 
response at the end of the combustor where only a few drops remain to burn and 
the combustion response is very dependent on when in the pressure cycle these 
remaining drops are vaporized. 

The gas oscillations and combustion response field repeat around the chamber 
with each baffle compartment. Therefore, only the first four tangential 
locations, corresponding to one baffle compartment, are shown in Fig 2 and 
discussed herein. Profiles inside the hub baffle (locations 2,2 4,1 and 4,4) 

are presented in Figure 2 as well as within a baffle compartment (locations 5,1 
5,4 7,2 8,1 and 8,4). 

The response for radial positions 1 through 3 (near the center and within the 
hub) and all tangential positions are very similar as the oscillations are 
spinning. These positions are represented by the response at location 2,2. Next 
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to the hub (radial position 4) the influence of the radial baffles outside the 
hub are noticed and the responses corresponding to the cells near the radial 
baffles (4,1 and 4,4) are shown. Within the baffle compartment the responses 
are given for the corners (5,1 5,4 8,1 and 8,4) and the center (7,2) to show 

the extreme variations noticed. 

At the first axial position, the real part of the response is always a large 
positive number and the imaginary part is a larger negative number. This is a 
result of the large variation in drop size produced by the oscillation in fuel 
flow during a cycle. In the baffle compartments, the drop size varied between 
270 and 950 microns during a cycle. With the large negative fuel response listed 
above the axial velocity producing atomization is very low when the pressure is 
high. Thus large drops are being produced when the pressure is high, resulting 
in low burning rates for these drops. Within the first cell 40 to 50 % of the 
oxidizer is burning (depending on the oscillation amplitudes) . The average time 
for the oxidizer to travel from the injector face (atomization plane) to the 
center of the first axial cell is 0.40 of a cycle. Therefore, "on average" the 
large drops produced at the high pressure portion of the oscillations are 
burning at the low pressure portion of the oscillation in the first cell. This 
produces the large "positive" combustion response in the first cell. The fact 
that the combustion response has a larger negative imaginary response indicates 
that the time delay has resulted in the large drops burning closer to the period 
when the pressure is near zero and increasing. 

The real part of the combustion response in the first cell varies between 1.78 
(2,2 position) to 3.41 (5,4 position). The imaginary part of the combustion 
response varies between -5.30 (8,2 position) and -7.67 (4,1 position). These 
combustion responses could not be considered a constant as assumed in the 
simple combustion models discussed above. 

The normalized combustion response drops very fast after the first cell, due to 
the decrease in burning rate axially and to the fact that the drops take 
different amounts of time to arrive at downstream cells. The effects of the 
variation in drop size with time are being "smeared" by the drops having 
different life times. It takes several cycles for the drops to arrive at these 
downstream cells. As a result, the contribution of the oscillating burning 
downstream is minimized. However, the variation in combustion response with 
radial and tangential directions are very large. Again, this violates the 
assumption used in the simple combustion models that the combustion response is 
uniform throughout the chamber. 

To provide more detail on the variation of combustion response with radial and 
tangential positions, a normalized stream tube response was determined. The 
normalized stream tube response was calculated by summing the local combustion 
rates. The variation in normalized stream tube response is given in Figure 3. 

Figure 3. clearly shows the large variation in combustion response with 
tangential position in the outer baffle compartments and with radius inside the 
hub compartment. The five blades of the baffle are symmetrically located around 
the chamber (one is at the 3 o'clock position). The combustion response is 
highest in the cells downstream of the baffle (at 2 o'clock) and is lowest in 
the cells blocking the baffle (at 12 o'clock). The highest combustion responses 
occur inside the hub. The phase angle of the response is indicated by a vector 
in Fig. 3 (zero phase angle occurs at 3 o'clock). The phase angle of the 
response is relatively constant in the bladed baffle compartments. A 45° phase 
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shift occurs between the hub compartment and the outer bladed compartments. 
Although the magnitude of the combustion response is highest in the hub 
compartment, the "poor" phase angle results in less response that is in phase 
with the pressure oscillation. The combustion zones in the bladed compartments 
are more important in driving oscillations in the SSME. 

To determine the best representation of a combustion response for the entire 
chamber the burning rates were summed over the full chamber volume. For the 
baffled chamber, the "normalized full chamber combustion response" was 1.358, - 
2.084i. This is in contrast to the value of 1.35, +0.0i determined using the 
FTE approach to drive the IT mode. Obviously, the combustion response will not 
drive the oscillations at the "natural frequency" since the phasing of the two 
vectors are not the same (the tuning frequency would be lower as indicated by 
Fig 1 for negative imaginary response). In addition, the local combustion 
oscillations calculated in the two stages are very different which means we do 
not have the correct oscillation profiles. 

Similar local combustion response calculations were made for the unbaffled SSME 
with no absorber. For the IT mode in an unbaffled chamber, all the angular 
positions see the same oscillation. Therefore, all the combustion responses are 
the same for a given radial position. The normalized stream tube combustion 
responses are shown in Table IV for the SSME with stability aids removed. 

TABLE IV 

STREAM TUBE COMBUSTION RESPONSE "Unbaffled Chamber" 


Response Ax Pos *1 


Rad 

Real 

Imaginary 

Real 

Imaginary 

1 

0.784 

-2.574 

1.673 

-3.951 

2 

0.294 

-3.235 

0.866 

-5.606 

3 

0.787 

-3.119 

1.682 

-5.807 

4 

1.139 

-2.909 

2.240 

-5.522 

5 

1.302 

-2.655 

2.567 

-5.045 

6 

1.336 

-2.455 

2.736 

- 4.738 

7 

1.329 

-2.303 

2.807 

-4.593 

8 

1.346 

-2.192 

2.891 

-4.575 


Normalized Full Chamber Combustion Response = 1.272, -2.453i 

Table IV again shows the large variation in combustion response with position 
with an unbaffled chamber. Since the wave is spinning, there are variations in 
only the radial and axial directions. The first axial position again has a very 
large value (twice the normalized stream tube value) and then rapidly drops, as 
seen in Fig. 2 for the baffled chamber. The largest response values are at the 
outside of the chamber where they are most effective in driving an oscillation. 
The calculated response for the full chamber has a slightly lower real value for 
the unbaffled chamber compared to the baffled chamber. Normally, this would 
indicate that the unbaffled chamber is "more stable" than the baffled chamber. 
With the SSME, the large combustion response is due to the varying atomization 
produced by the large injector fuel flow response. 

The results and discussions above indicated that the fuel response is 
responsible for a large portion of the combustion response in the SSME. This 
result was checked by performing calculations with different portions of the 
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combustion model turned off. Calculations were also performed with different 
oscillations amplitude levels to show the nonlinear characteristics of the 
combustion process. The results of these calculations are shown in Table V. 

TABLE V 

FACTORS INFLUENCING COMBUSTION RESPONSE 


Combustion Response 


Factor 

5% Pk- 

Pk Amp 

20% 

Pk-Pk 

80% 

Pk-Pk 


Real 

Imag 

Real 

Imag 

Real 

Imag 

Fuel Flow Osc 

1.453 

-3.505 

0.803 

-2.268 

0.166 

-0.909 

Axial Vel. Osc 

0.025 

0.652 

0.025 

0.655 

0.020 

0.424 

Den. Osc 

0.419 

-0.046 

0.417 

-0.069 

0.405 

-0.085 

Ox Flow Osc 

0.242 

-0.322 

0.277 

-0.341 

0.417 

-0.232 

Pres. Osc 

0.070 

0.010 

0.064 

0.011 

0.045 

0.000 

Vr & Va Osc 

0.000 

0.000 

0.001 

-0.001 

-0.008 

-0.002 

Sum of All 

(2 . 209) ( 

-3.201) 

(1.414) (-2.044) 

(1 . 045) ( -0 . 804) 

All Osc 

2.223 

-3.026 

1.358 

-2.084 

0.646 

-1.287 


Table V confirms the importance of fuel flow oscillations on the combustion 
response. It is the phenomena making the largest contribution to combustion 
response at all amplitudes. The influence is very nonlinear, decreasing with 
increasing amplitude. At higher oscillation amplitudes the velocity of the fuel 
in the injector element reaches sonic velocity during a portion of the 
oscillation. Since the velocity cannot exceed sonic velocity, this reduces the 
fuel flow oscillation. To account for this fact, the fuel response in HICCIP is 
decreased by -the square root of the oscillation amplitude at which sonic is 
first achieved divided by the actual amplitude. For the 20% pk-pk chamber 
pressure oscillations , sonic is reached in the four rows within the baffle 
compartment and the response is reduced to 0.8 of its value at low amplitude. At 
80% Pc pk-pk, rows 2 through 10 reach sonic and the fuel response in the baffle 
compartment area is reduced to 0.4 of its value at low amplitude (values in 
Table III). As a result of this feature, the contribution of the fuel flow 
response decreases with increasing amplitude. The other terms are relatively 
insensitive to the oscillation amplitude at low amplitudes but do change at very 
high amplitudes. 

Axial velocity oscillations are the next most important influence on combustion 
response after the fuel response. With the presence of the baffles, an axial 
wave is present in the baffle compartments, thereby producing large axial 
velocity oscillations. With a steady axial velocity in all cells, the axial 
velocity produces a sinusoidal oscillation in burning rate ( the burning is 
Reynolds Number dependent) . The baffles also reduce the frequency of the 
oscillation. Within the hub, the reduced frequency combined with the smaller 
diameter results in a large amplitude increase as you move from the injector to 
the tips of the baffle. This is accompanied by large axial velocity oscillations 
within the hub. The contribution of the axial velocity is out of phase 
(imaginary response) as the axial velocity oscillations are generally 120 
degrees out of phase with the pressure. 

Density oscillations always produce a 0.4 contribution to the combustion 
response. This is associated with the density term in the Reynolds number which 
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influences the burning rate in proportion to the square root of the Nusselt 
number. Since density is proportional to pressure divided by gamma, the 
predicted influence of density on response is about 0.4 and is independent of 
amplitude, frequency, drops size, etc. 

Oxidizer flow oscillations have a nominal influence on the combustion response 
in the SSME. Oxidizer injector pressure drop is not excessive and as Table II 
shows the oxidizer flow response is as large as -1.2 at some rows. The negative 
flow response results in a positive combustion response because the time lag 
between injecting the oxidizer and burning is close to a half cycle. The 
oxidizer flow response is reduced due to the "smearing" of the flow oscillations 
through many cycles in the combustion process. The influence of the oxidizer 
flow oscillations increases with amplitude due to the change in drop size and 
velocity, which changes the time lag between injection and burning as well as 
the burning rates . 

Radial and tangential gas velocity and pressure oscillations have little 
influence on the combustion response. Without steady state or time averaged 
radial and tangential velocities, the effect of these oscillations on burning 
rate cancel (burning rate is a function of the magnitude of velocity) . The 
radial and tangential gas velocities would contribute to the harmonic terms in 
the Fourier expansion. The pressure oscillation also has little influence as it 
changes the "surface" temperature of the drop which slightly changes the heat 
required to vaporize the drop. Pressure also slightly influenced the drop size 
and atomization rate via the density influence in the atomization process which 
was considered as a pressure oscillation influence. These factors were all small 
and not very dependent on oscillation amplitude. 

Contributions of the various factors toward the total combustion response is 
additive at the lower amplitudes. The term "Sum of All" in Table V is the value 
obtained by adding all the individual contributions. At the 5 and 20 % Pk-Pk 
levels, the actual calculated response using all the terms is within 5% of the 
individual sums . At higher amplitudes, the factors interact, as one would 
expect, and the sum appears to be quite different than the actual. The major 
influence however has been to change the phasing of the combustion response 
(larger imaginary and lower real values actually). Basically, this means that 
the time lag has changed. 

All of the results presented above are very dependent on the constants used in 
the atomization process to determine drop size and atomization rates. These 
will change the steady state or time averaged profiles as well as the time lags 
between injection and burning. To show this influence, calculations were made 
with the same gas phase oscillations but with different atomization constants. 
These results are shown in Table VI. 
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TABLE VI 


INFLUENCE OF ATOMIZATION CONSTANT ON COMBUSTION RESPONSE 


Combustion Response 


CONSTANTS 

DrSz 

50% 

Baffled 

Unbaffled 

Size 

Rate 

Micron 

Vap 

Real 

Imag 

Real 

Imag 

18800 

0.2666 

268 

0.4 

0.46 

-4.70 

0.21 

-4.74 

26320 

ft 

375 

0.5 

1.11 

-3.29 

0.96 

-3.60 

37600 

tt 

507 

0.9 

1.36 

-2.08 

1.27 

-2.45 

52640 

ft 

711 

1.6 

1.43 

-1.00 

1.40 

-1.59 

75200 

tf 

1072 

2.8 

1.21 

-0.06 

1.18 

-0.79 

37600 

0.1333 

1558 

5.7 

2.59 

-1.30 

2.63 

-1.30 

rt 

0.1866 

712 

1.5 

1.52 

-2.15 

1.36 

-2.83 

n 

0.2666 

507 

0.9 

1.36 

-2.08 

0.96 

-2.45 

Tf 

0.5333 

480 

0.8 

1.39 

-1.84 

1.43 

-2.07 


Table VI shows the strong "tuning" influence of the atomization constants. Both 
the size and rate constants show "optimum” values at which the response is 
either a maximum or minimum. The size constant of 37,600 (used in all the 
calculations presented above) results in a response for both the baffled and 
unbaffled chamber that is almost a maximum. Increasing or decreasing the size 
would improve stability. The rate constant of 0.2666 ( also used above) results 
in a "minimum" response for both the baffled and unbaffled chamber. Drop size 
and location of the zone where 50% of the oxidizer is vaporized is also shown. 
Changing the rate constant changes the atomization plane and the drop size. The 
change in drop size is due to a change in the gas velocity field surrounding the 
liquid when the atomization rate is changed. 

Table VI also shows the importance of the atomization constants in comparing the 
stability of the baffled vs unbaffled chamber. With either larger or smaller 
rate constants the baffled chamber is more stable than the unbaffled. At the 
rate constant used in the previous results, the unbaffled chamber is more stable 
than the unbaffled. At smaller size constants, the unbaffled chamber becomes 
much more stable than the baffled. 


COMPUTATION TIMES 

The gas oscillation and HICCIP vaporization models have been used to determine 
the computer time required for different grid sizes. An example of the 
computation time using a PC with an 80486 (33MHz) cpu to determine a solution 
for the gas oscillations with a specified oscillating burning rate and the 
combustion response with a specified gas oscillation is shown in Table VII. 
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TABLE VII 

COMPUTATION TIME FOR SINGLE SOLUTION - 1 T MODE 
486 Computer 33Mhz 


Cell Size 
Ni Nj Nk 
1 1 16 

4 18 

8 1 16 

16 1 32 

4 10 8 

8 20 16 

8 20 32 


Model 


Computer Time 


Gas Oscillations 


Comb Response 
0.04 Min 
0.15 " 

0.3 

1.2 

3 

12 

12 


Baffled 


Cylinder 


II 


It 


It 


0.05 Min 
0.07 " 

2.30 " 

18.19 " 

0.03 Hrs 
1.62 Hrs 
8.03 Hrs 


These results indicate that the computer time for a reasonable representation of 
a baffled chamber (8 radial cells, 20 tangential and 32 axial) is not excessive. 
The time to calculate the combustion response is much less than required to 
determine the gas oscillations . For a simple 1-D chamber (1x1x16), the gas 
oscillations require the same amount of time as the combustion response 
calculations. The time for gas oscillations varies with the 1.5 power of grid 
size while the combustion response time varies directly with the radial and 
tangential grid size. The effective axial grid size for the combustion response 
is always much finer than used in the gas oscillations program and is determined 
by the frequency and distance required to vaporize each drop size. 

After a computation to determine the flow oscillations with a specified 
oscillating burning rate and a burning rate oscillation with that flow 
oscillation you find that the burning rate is not consistent with oscillations 
as one would expect from physical models. Therefore, the process has to be 
repeated several times (3 to 10 depending on the initial guess) before the 
burning rate agrees with the oscillations. Typical computational times for 
completing this type of calculation to determine the combustion response at a 
given frequency (equivalent to determining the complex frequency with a given 
combustion model) is: 


Computation time is dependent on the initial guess for the burning rate 
oscillations. If a very poor guess is used many more iterations are required 
before the burning rate oscillations agree with the calculated pressure profile 
and constant combustion response assumption. If a frequency scan is used and the 
solution from the last frequency is used as the next assumption, the above 
computer times can be decreased by as much as a factor of 4 (four) . Similarly 
higher order modes and/or frequencies require significantly less computation 
time (as much as 1/3 for the 3T mode). Computer times are smaller when the 
original guess for the pressure and velocity perturbations in the Conjugate 


TABLE VIII 


COMPUTATION TIME FOR RESPONSE (FREQUENCY) SOLUTION 
Baffled Chamber IT Mode 


Cell Size 
Ni Nj Nk 
4 10 8 

8 20 16 

8 20 32 


COMPUTER TYPE 

486 (33Mhz) VAX Cray 

6 . 6 Min 1. 2 Min 

3.19 Hrs 36 Min 1.55 Min 

13.7 Hrs 1.3 Hrs 3.47 Min 
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Gradient Solver are zero. A good initial guess for the burning rate oscillation 
is important. 


SUMMARY 

A brief review of instability modeling prior to 1972 was made. Current 
approaches to instability modeling were outlined and the advantages and 
disadvantages of each approach were discussed. The theory for the Fourier time 
expansion approach was outlined in detail. Calculations were made for the 
combustion response profiles in the SSME combustor (with a five bladed hub 
baffle and absorbers) . The contributions of various physical phenomena 
(atomization, injection, etc.) to the combustion response profiles were assessed. 
Computational requirements for the Fourier time expansion were discussed. 
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SSME Response Curves 
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Figure 2. - Normalized local combustion responses at 

various radial and tangential grid positions 
( 8 radial x 20 tangential x 30 axial cells ) 
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Figure 3. - Normalized stream tube responses for a baffled chamber 
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